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(A)  Windowed  waveforms  (1  second  before  and  2.82  seconds  after  the  catalog  P  pick) 
for  two  New  Zealand  earthquakes  recorded  at  station  MOW.  The  station  sampling 
rate  is  50  sample/sec  and  the  signals  are  aligned  by  the  P  picks.  (B)  Windowed 
waveforms  band-pass  filtered  between  1  and  6  Hz.  (C,  D)  CC  series  for  the  (C) 
unfiltered  (raw)  and  (D)  filtered  waveforms.  (E,F)  Bispectrum-correlation  series  for 
the  (E)  raw  and  (F)  filtered  waveforms.  (G)  Stacked  raw  waveforms.  The  upper  trace 
is  obtained  by  shifting  the  waveform  of  the  first  event  with  the  CC  determined  lag 
relative  to  the  second  one,  while  the  lower  trace  is  computed  after  shifting  the 
waveform  of  the  first  event  with  the  lag  calculated  with  the  bispectrum  method.  The 
root-mean-square  amplitudes  for  the  two  traces  are  0.25  and  0.27  respectively.  (H) 
Same  as  G  for  the  stacked  band-pass  filtered  waveforms.  The  root-mean-square 
amplitudes  for  the  two  traces  are  0.32  and  0.34  respectively. 
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Figure  2. 


Examples  of  a)  rejection  of  a  high  CC  value  due  to  inconsistent  results  between  the 
CC  and  BS  lags  and  b)  acceptance  of  a  low  CC  value  due  to  consistent  results  between 
the  CC  and  BS  lags. 
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Figure  3. 


Locations  of  53  earthquakes  near  lake  Wairarapa,  NZ.  (A)  Before  relocation;  (B) 
relocated  with  catalog  differential  travel  times;  (C)  relocated  with  CC  differential 
travel  times  chosen  with  the  threshold  criterion;  (D)  relocated  with  CC  differential 
travel  times  verified  with  the  bispectrum  method;  (E)  relocated  with  both  catalog  and 
threshold-chosen  CC  differential  travel  times;  (F)  relocated  with  both  catalog  and 
bispectrum-verified  CC  differential  travel  times. 

Figure  4.  23 

Maps  of  (a)  reference  events  and  (b)  test  events  and  test  stations  from  our  California 
ground-truth  dataset. 

Figure  5.  23 

Contour  plots  (map  view  and  depth  sections)  of  the  maximum  cross-correlation  value 
per  bin  for  the  7  stations  for  a  test  event  (located  at  lat.  34.3°,  Ion.  -1 18.58°,  and 
depth  1.8  km,  cross)  correlated  against  the  set  of  reference  events  in  the  (a)  low- 
frequency  band  (step  one)  and  (b)  high-frequency  band  (step  two).  Note  the  good 
low-resolution  location  estimate  on  the  left  (circle)  and  the  improved  high-resolution 
location  estimate  (circle)  on  the  right. 
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Figure  6. 

Example  application  of  KLD  to  a  set  of  seismograms  with  a  pair  of  "arrivals"  only 
one  of  which  is  time-aligned  across  the  records  (top  left),  with  noise  added  (top  right). 

The  data  reconstructed  from  the  6 1  eigenvectors  with  eigenvalues  above  average  in 
size  (bottom  left)  retains  both  the  time-aligned  arrival  and  the  arrival  with  moveout 
(middle  left),  leaving  little  signal  energy  in  the  residual  image  (middle  right).  In 
contrast,  eigenimage  decomposition  cannot  reconstruct  both  arrivals  adequately. 

Figure  7.  25 

Example  application  of  KLD  to  real  data  (upper  left).  Large  and  small  eigenvalues 
were  excluded  in  this  case.  The  reconstructed  data  (upper  right)  from  the  30 
eigenvectors  with  eigenvalues  from  #5  to  #34  (middle  right)  produces  similar  looking 
waveforms,  leaving  relatively  little  energy  in  the  residual  image  (middle  left).  The 
eigenvectors  for  the  larger  eigenvalues  (lower  left)  and  their  corresponding 
reconstructions  (lower  right)  show  the  long-period  noise  character  of  the  first  4 
eigenvectors. 

Figure  8.  26 

a,b)  Seismic  data  windowed  around  the  P  arrival  for  the  a)  shallow  (2  km)  and  deep 
(12  km)  earthquakes  from  Northridge  (S.  California)  cluster.  The  shallow  event  is  used 
as  the  reference  signal  and  the  deeper  event  is  used  as  the  primary  signal;  noise  has 
been  added  to  the  latter  trace,  c)  Zero-lag  cross-correlation  between  the  primary  and 
CANC-filtered  traces  for  a  range  of  signal  alignments  (see  text),  d)  SNR  for  the 
CANC-filtered  trace  versus  the  error  trace  for  a  range  of  signal  alignments,  e)  Common 
signal  obtained  from  CANC  adaptive  filter,  f)  Error  trace.  Coherency  filtering  applied 
to  data  in  a)  and  b)  results  in  the  signal  traces  displayed  in  g)  and  h),  respectively. 

Figure  9.  27 

A  horizontal  slice  through  the  true  synthetic  velocity  model.  The  true  velocity  model 
in  3D  is  similar  to  a  "vertical  sandwich,"  with  velocity  constant  with  depth. 

Figure  10.  27 

Event  locations  (filled  circles)  and  stations  (triangles)  used  for  the  synthetic  data  set. 

The  inversion  grid  used  in  the  standard  and  DD  tomography  solutions  is  shown  as 
crosses.  The  inversion  grid  points  are  at  X=-35,  -15,  0,  2,  4,  6,  20,  35  km,  at  Y=-60,  - 
40  -20,  0,  20,  40  km,  and  at  Z=0,  3,  7,  11,  16  km. 

Figure  11.  27 

Horizontal  slices  through  the  velocity  models  from  (a)  standard  tomography  and  (b) 

DD  tomography,  and  the  difference  between  (c)  the  DD  tomography  solution  and  the 
true  model  and  (d)  the  standard  tomography  solution  and  the  true  model.  Black  dots 
indicate  the  earthquake  hypocenters  within  half  the  grid-size  of  the  slice. 
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Figure  12. 

Comparison  of  earthquake  locations  along  the  Hayward  fault  determined  using  (a)  DD 
tomography,  (b)  DD  location,  (c)  standard  tomography,  and  (d)  catalog  data  (no  DD 
data).  Each  part  shows  a  lat-lon  plot  (upper  left)  and  zoomed-in  plots  (rotated  to 
fault-parallel/fault-normal)  of  epicenters  (upper  right)  and  fault-normal  (lower  left) 
and  fault-parallel  (lower  right)  depth  sections.  Note  how  closely  the  hypocenters  in 
panels  a  and  b  resemble  each  other,  but  with  a  systematic  SW  shift.  The  tight 
clustering  is  not  evident  in  either  the  standard  tomography  solution  (c)  or  the  catalog 
locations  using  absolute  picks  only  (d). 

Figure  13.  29 

(a)  The  left-hand  panels  show  fault-normal  slices  through  the  DD  tomography  model. 

The  model  has  sharper  features  and  a  better  correspondence  between  the  main  region 
of  seismicity  and  the  position  of  the  sharp  velocity  contrast  compared  to  the  standard 
results  in  b.  (b)  The  right-hand  panels  show  fault-normal  slices  through  the  standard 
tomography  model.  The  model  has  features  that  are  less  sharp  and  a  poorer 
correspondence  between  the  main  region  of  seismicity  and  the  position  of  the  velocity 
contrast  than  the  DD  results 

Table  1.  30 

The  absolute  differences  between  the  true  locations  and  those  from  the  DD  location 
method  based  on  1 D  velocity  model,  standard  tomography,  and  DD  tomography. 
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Summary 


Our  research  program  consists  of  four  components,  each  involving  some  aspect  of 
multiple-event  analysis:  (1)  high-precision  waveform  cross-correlation  (WCC)  for  arrival  time 
estimation,  (2)  robust  event  clustering,  (3)  waveform  decomposition  and  source  wavelet 
deconvolution  reshaping,  (4)  double-difference  (DD)  multiple-event  location  and  tomography. 
Our  research  focused  initially  on  the  development  and  testing  of  these  seismic  analysis  methods 
using  "ground-truth"  datasets  at  different  scales  (local  and  regional),  and  was  followed  by  the 
application  of  these  methods  to  real-time  and  simulated  real-time  data  streams. 

We  completed  the  development  and  carried  out  initial  testing  of  a  multiple  eigentaper 
algorithm  for  improved  integer-sample  WCC.  This  is  an  extension  of  the  use  of  the  eigentaper 
approach  from  the  sub-sample  to  the  integer-sample  computations.  Modifications  have  been 
made  to  the  WCC  code  logic  to  provide  the  initial  step  for  adaptation  for  multiple  frequency- 
scale  correlations.  In  addition,  modifications  have  been  completed  to  interface  the  WCC  code 
with  our  new  DD  tomography  code.  This  alters  the  code  output  to  produce  station  correlations 
and  lags  on  a  by-event  basis.  We  initiated  an  investigation  of  the  use  of  multi-wavelet 
decomposition  methods  (in  place  of  Fourier  decomposition)  for  potential  improvement  of  WCC. 
The  initial  findings  were  promising,  but  further  assessment  is  required. 

Differential  arrival  times  for  pairs  of  seismic  events  observed  at  the  same  station  are  often 
calculated  by  WCC  techniques.  Researchers  generally  choose  the  differential  times  to  use  based 
on  the  associated  cross-correlation  (CC)  values  exceeding  a  specified  threshold.  When  two  similar 
time  series  are  corrupted  by  correlated  noise,  the  time  delay  estimate  calculated  with  the  WCC 
technique  may  not  be  reliable.  Thus  the  selection  of  a  threshold  value  is  important.  If  it  is  set  too 
high,  then  only  a  limited  number  of  very  accurate  differential  time  data  are  available  to  constrain 
the  relative  positions  of  earthquakes.  If  the  threshold  value  is  set  too  low,  then  many  unreliable 
differential  time  estimates  are  used  and  they  will  negatively  affect  the  final  relocation  results.  The 
bispectrum  method  can  suppress  the  correlated  Gaussian  noise  sources  in  two  similar  time  series 
and  be  used  to  obtain  the  relative  time  delay  between  them.  We  compute  time  delay  estimates 
between  the  waveform  pairs  with  the  bispectrum  method  and  use  them  to  verify  the  WCC- 
determined  one.  This  technique  can  provide  quality  control  over  the  selected  time  delay  estimates 
and  potentially  provide  more  differential  travel  time  data  for  close  events  pairs,  by  verifying  the 
reliability  of  differential  times  that  might  not  meet  the  threshold  criterion.  A  manuscript 
describing  the  technique  and  illustrating  its  application  has  been  published  (Du  et  al.,  2004).  We 
submitted  the  CC-bispectrum  analysis  code  (BCSEIS)  to  our  Product  Integrator. 

We  have  tested  two  versions  of  a  two-step  event  location  technique  with  California  and 
New  Mexico  ground  truth  (GT)  data  set.  In  step  I  of  the  first  (UW)  version,  event  locations  are 
estimated  via  WCC  of  Hilbert-envelope  waveforms.  The  Step  1  procedure  accurately  located  the 
epicenters  within  20  km  (approximately  1  bin)  of  the  NCSN  catalog  location  for  16  of  the  20  test 
events.  However,  depth  estimates  were  poorly  constrained.  The  Step  2  refinement  based  on 
WCC  of  the  high-frequency  P  arrivals  followed  by  DD  relocation  improved  both  epicenter  and 
depth  estimates.  This  two-step  procedure  offers  advantages  over  conventional  location 
techniques  including  a  Step  1  event  location  not  dependent  on  the  accuracy  of  catalog  picks  and  a 
Step  2  improvement  in  location  accuracy  using  the  high-accuracy  DD  technique.  For  the  other 
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(NMT)  version,  the  approach  is  to  find  the  most  similar  event  from  a  reference  database,  and 
assign  phases  and  location  information  of  the  most  similar  event  in  a  reviewed  database  according 
to  the  spectrogram  cross-correlation.  No  preliminary  pick  or  location  information  is  needed  for 
the  unknown  event.  Both  approaches  have  proven  effective  in  our  trials. 

We  developed  a  wavelet-based  auto-picker  for  detecting  and  picking  first-P  arrivals  that 
will  be  used  in  conjunction  with  our  location  processing  procedures,  both  for  identifying 
waveform  segments  for  correlation  and  for  adjusting  stack  picks.  A  manuscript  describing  the 
approach  has  been  published  (Zhang  et  al.,  2003). 

We  tested  the  potential  of  waveform  reshaping  methods  using  the  California  GT  dataset. 
Two  approaches  were  examined:  an  empirical  deconvolution  method  and  a  deterministic  depth- 
based  method.  The  key  to  both  methods  is  to  determine  an  approach  that  "simplifies"  the  first- 
arrival  wavelet  but  preserves  the  relative  arrival  time. 

We  tested  a  variety  of  decomposition  and  filtering  techniques  with  the  multiple  goals  of 
noise  reduction,  waveform  "homogenization"  (reshaping  to  similar  "wavelet"),  and  depth  phase 
identification.  Methods  examined  include  eigenimage  decomposition,  Karhunen-Loeve 
decomposition  (KLD),  Correlated  data  Adaptive  Noise  Canceling  (CANC),  and  coherency 
filtering.  A  key  issue  is  the  effect  of  arrival  misalignment  between  a  pair  of  traces  on  the  success 
of  a  given  technique.  We  tested  a  two-stage  process  whereby  the  similarity  and  strength  of  a 
CANC-filtered  trace  (relative  to  the  original  and  error  traces,  respectively)  is  used  to  align  the 
trace  to  the  reference  signal. 

We  completed  the  development  and  testing  of  a  3D  DD  tomography  code,  starting  from 
the  DD  location  code  hypoDD.  We  applied  our  new  DD  tomography  algorithm,  tomoDD,  to  a 
dataset  from  the  Hayward  fault,  California,  as  an  initial  test  to  compare  to  previously  published 
location  results.  We  published  a  manuscript  on  the  method  and  initial  application  (Zhang  and 
Thurber,  2003).  We  also  submitted  the  DD  tomography  code  (tomoDD)  to  our  Product 
Integrator.  An  additional  version  of  the  DD  tomography  code  was  subsequently  developed,  a 
spherical-earth  version  appropriate  for  regional-scale  applications. 
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Introduction 


The  standard  processing  approach  for  the  determination  of  seismic  event  hypocentral 
locations  involves  a  single-event  procedure.  Phase  arrivals  are  timed  and  associated  to  define  an 
event,  and  the  associated  phases  are  used  (potentially  along  with  arrival  azimuths  and 
slownesses)  to  compute  a  best-fitting  location.  It  has  been  shown  that  order-of-magnitude 
reductions  in  event  location  errors  and  uncertainties  are  possible  if  multiple-event  techniques  are 
utilized.  Recent  studies  of  earthquakes  in  Hawaii  (Rubin  et  al.,  1998),  California  (Waldhauser  et 
al.,  1999),  and  at  the  Soultz  geothermal  area  (Rowe  et  al.,  2002b),  and  of  explosions  at  the 
Balapan  test  site  (Phillips  et  al.,  2001;  Thurber  et  al.,  2001),  and  others,  have  demonstrated  in  a 
dramatic  fashion  the  substantial  improvement  in  the  definition  of  seismogenic  features  or  in  the 
accuracy  of  relative  locations  of  ground-truth  events  that  is  possible  with  multiple-event  methods 
utilizing  high-precision  arrival-time  estimation. 

The  accuracy  of  event  hypocenters  is  determined  by  several  factors,  including  the 
network  geometry,  available  phases,  and  arrival  time  accuracies  (Pavlis,  1986).  If  more  phases  are 
available,  the  event  hypocenters  may  be  more  stable.  For  example,  S-wave  phases  will  be  very 
useful  in  constraining  event  depths  and  providing  information  that  helps  to  decouple  the 
hypocenters  from  the  structure  in  the  inversion  (Gomberg  et  al.,  1990).  Due  to  the  presence  of 
noise,  the  arrival  times  picked  either  manually  or  automatically  generally  have  errors.  Recent 
studies  have  shown  spectacular  improvements  in  location  precision  for  earthquakes  and 
explosions  when  WCC  and  event  clustering  techniques  are  used  to  improve  arrival  time  estimates 
or  determine  high-precision  relative  arrival  times  (VanDecar  and  Crosson,  1990).  Such  studies  are 
based  on  the  assumption  that  waves  generated  by  two  similar  sources,  propagating  along  similar 
paths,  will  generate  similar  waveforms,  and  WCC  can  then  be  used  to  determine  precise  relative 
arrival  times.  These  studies  have  demonstrated  substantial  improvement  in  the  definition  of 
seismogenic  features  and  in  the  accuracy  of  relative  locations  of  ground-truth  events  that  is 
possible  using  multiple-event  methods  with  high-precision  absolute  or  relative  arrival-time  data. 

Most  seismic  stations  exhibit  signal  quality  problems  to  some  extent  or  another.  When 
analyzing  datasets  that  cover  a  wide  geographic  expanse  and  a  wide  magnitude  range,  any  master 
station  may  at  times  have  poorer  signal  quality  than  other  sensors  in  a  network,  either  because  of 
clipping  for  a  large,  nearby  event,  or  poor  signal/noise  for  smaller,  more  distant  ones.  Reliance  on 
a  single  station  to  perform  all  similarity  choices  is  therefore  not  generally  the  best  approach,  and 
we  are  testing  the  success  of  a  more  adaptive  technique  where  two  or  more  stations  may  have 
their  waveform  similarity  matrices  compared  for  consistency,  particularly  among  events  that 
have  been  orphaned  in  the  clustering  step  for  one  of  the  stations. 

Often  the  first  arrival  is  very  emergent  and  only  later  portions  of  the  signal  can  be 
compared.  In  these  cases,  it  is  useful  to  have  a  sliding  correlation  window  to  search  for  similarity 
later  in  the  event  coda.  We  have  performed  interactive  testing  to  demonstrate  that  sometimes 
highly  similar  events,  whose  first  arrivals  were  uncorrelatable,  can  be  well  aligned  using  the 
secondary  phases.  One  does  not  generally  expect  even  among  similar  events  that  alignment  of 
later  phases  will  necessarily  optimally  align  the  first  arrivals;  however,  among  events  whose 
similarity  has  been  established  at  several  correlation  lengths  for  several  stations,  this  provides  a 
means  of  identifying  arrivals  at  noisier  stations  and  improving  the  ability  to  estimate  hypocenters 
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for  these  events.  We  see  this  approach  as  a  potential  tool  in  generating  autopicks  at  stations  with 
poorer  S/N;  if  a  new  event  has  been  identified  as  belonging  to  a  family  of  highly  similar  events, 
the  pick  for  the  poor  station  may  be  predicted  based  on  its  picks  for  better-recorded  events  in  the 
same  family,  then  the  validity  of  the  autopick  may  be  verified  (and  its  time  adjusted)  through 
comparing  later,  stronger  portions  of  the  waveforms. 

WCC  analysis  of  P-wave  data  has  shown  the  ability  to  identify  event  clusters  based  on 
the  cross-correlation  information.  In  a  small  study  region,  P  wave  data  might  suffice  for 
identifying  clusters  and  hence  the  source  position.  However,  the  smaller  window  lengths  and 
high-frequency  data  used  in  WCC  P-wave  analysis  limits  the  ability  to  construct  unique  location 
information  in  a  larger  study  region.  Shearer  ( 1 994)  has  demonstrated  the  use  of  low-frequency 
data  (3  mHz  to  0  0.1  Hz)  for  global  event  location  and  Withers  et  al.  (1999)  have  successfully 
applied  a  similar  method  to  local  and  regional  data  utilizing  coarse  and  fine  spatial  grids.  A 
hierarchical  clustering  analysis  incorporating  both  the  statistical  advantages  of  cluster  analysis 
with  adaptive  use  of  grid  size  and  frequency  content  of  the  data  set  offers  an  alternate  means  to 
locate  earthquakes.  The  order  of  magnitude  reductions  in  location  errors  found  with  the  use  of 
WCC  analysis  for  event  clusters  suggest  a  similar  improvement  can  result  from  our  hierarchical 
clustering  analysis,  with  location  accuracy  limited  mainly  by  the  distribution  of  the  ground  truth 
(reference  event)  locations,  binning  interval,  and  the  clustering  statistics. 

Differential  arrival  times  for  pairs  of  seismic  events  observed  at  the  same  station  are  often 
calculated  by  WCC  techniques.  These  differential  times  (or  adjusted  picks  derived  from  them) 
can  be  used  to  improve  dramatically  the  earthquake  relocation  results  (Got  et  al.,  1994;  Dodge  et 
al.,  1995;  Shearer,  1997;  Rubin  et  al.,  1999;  Waldhauser  and  Ellsworth,  2000;  Schaff  et  al.,  2002). 
Researchers  often  choose  the  differential  times  to  use  based  on  the  associated  CC  values 
exceeding  a  specified  threshold.  For  example,  Schaff  et  al.  (2002)  only  selected  those  time  delays 
with  CC  values  larger  than  0.7  and  mean  coherences  above  0.7.  The  selection  of  an  "optimum" 
threshold  value  is  important  but  difficult.  If  it  is  set  too  high,  then  only  a  limited  number  of  very 
accurate  differential  time  data  are  available  to  constrain  the  relative  positions  of  earthquakes.  On 
the  other  hand,  if  the  threshold  value  is  set  too  low,  then  many  unreliable  differential  time 
estimates  are  used  and  will  negatively  affect  the  final  relocation  results.  Also  the  time  delay 
estimate  calculated  with  the  CC  technique  may  not  be  reliable  when  the  noise  sources  in  the  two 
time  series  are  correlated. 

The  selection  of  an  optimum  threshold  value  is  important  but  difficult.  If  it  is  set  too 
high,  then  only  a  limited  number  of  very  accurate  differential  time  data  are  available  to  constrain 
the  relative  positions  of  earthquakes.  If  the  threshold  value  is  set  too  low,  then  many  unreliable 
differential  time  estimates  are  used  which  can  negatively  affect  the  relocation  results.  In  our 
approach,  we  compute  two  more  time  delay  estimates  between  the  waveform  pairs  with  the 
bispectrum  (BS)  method  (the  raw  data  and  band-pass  filtered  data)  and  use  them  to  verify  the 
WCC-determined  one  using  the  band-pass  filtered  data.  The  BS  method,  which  works  in  the 
third-order  spectral  domain,  can  suppress  correlated  Gaussian  or  low-skewness  noise  sources 
(Nikias  and  Raghuveer,  1987;  Nikias  and  Pan,  1988).  Du  et  al.  (2004)  adopt  this  method  to 
calculate  two  additional  time  delay  estimates  with  both  the  raw  (unfiltered)  and  band-pass 
filtered  waveforms,  and  use  them  to  verify  (select  or  reject)  the  one  computed  with  the  CC 
technique  using  the  filtered  waveforms.  This  technique  can  provide  quality  control  over  the 


4 


selected  time  delay  estimates  and  potentially  provide  more  differential  travel  time  data  for  close 
events  pairs,  by  verifying  the  reliability  of  differential  times  that  might  not  meet  the  threshold 
criterion. 

The  WCC  approach  and  its  variants  work  in  the  second-order  spectral  domain.  When  the 
underlying  signal  inside  two  similar  time  series  can  be  regarded  as  a  non-Gaussian  process  and  the 
noise  sources  as  zero-mean  Gaussian  processes,  the  similarities  between  the  two  time  series  can 
be  better  compared  in  the  third-order  or  bispectrum  domain  (Nikias  and  Raghuveer,  1987;  Nikias 
and  Pan,  1988;  Yung  and  Ikelle,  1997).  This  bispectrum  method  takes  advantage  of  the  fact  that 
for  Gaussian  processes  only  all  spectra  of  order  higher  than  two  are  identically  zero.  Thus  the 
effect  of  correlated  Gaussian  noise  is  completely  suppressed  when  we  compare  the  two  time 
series  in  the  bispectrum  domain,  while  it  may  make  WCC  techniques  fail  to  work  well.  Noise  at  a 
station  for  different  events  can  be  expected  to  be  partially  correlated  due  to  a  combination  of 
constant  predominant  noise  sources  with  time-varying  amplitude  (microseisms,  wind  or  cultural 
noise)  and  site  response  effects.  Cycle  skips  are  another  potential  WCC  pitfall  that  bispectrum 
analysis  may  help  detect. 

In  some  situations,  otherwise  similar  waveforms  may  be  "contaminated"  by  some 
interfering  signal  (e.g.,  a  depth  phase)  that  may  degrade  their  cross-correlation.  In  such 
circumstances,  it  would  be  desirable  to  be  able  to  extract  the  underlying  similar  waveforms, 
filtering  out  the  "noise"  of  the  contaminating  signal.  A  recent  paper  by  Ulrych  et  al.  (1999) 
reviews  several  decomposition  and  transformation  techniques  that  are  useful  for  signal-noise 
separation  in  an  exploration  seismic  context.  The  two  most  relevant  to  our  work  are  eigenimage 
decomposition  and  Karhunen-Loeve  decomposition  (KLD).  KLD  is  distinguished  from  the  more 
common  eigenimage  decomposition  (essentially  principal  components)  in  that  the  latter  extracts 
eigenvector-eigenvalue  information  from  multi-channel  data  based  on  a  covariance  matrix 
computed  at  a  single  lag  (zero),  whereas  the  former  carries  out  a  decomposition  for  single-channel 
data  based  on  a  covariance  matrix  for  multiple  lags  (Ulrych  et  al.,  1999).  In  practice,  to  apply 
KLD  to  multi-channel  data,  an  average  covariance  matrix  can  be  computed. 

We  can  make  an  important  distinction  between  the  two  fundamentally  different  ways 
WCC  data  has  been  used:  ( 1 )  by  directly  using  relative  arrival  times  to  determine  relative  event 
locations  (e.g..  Got  et  al.,  1994;  Waldhauser  and  Ellsworth,  2000),  and  (2)  by  adjusting  absolute 
arrival  time  picks  to  minimize  discrepancies  among  relative  arrival  times  (Dodge  et  al.,  1995; 
Rowe  et  al.,  2002a).  The  advantage  of  the  former  approach  is  that  it  incorporates  all  the  available 
information  contained  within  the  multitude  of  relative  arrival  time  differences  with  a  direct 
measure  of  quality  (the  correlation  value)  associated  explicitly  with  each  datum.  A  disadvantage 
is  that  some  simplifying  assumption  needs  to  be  made  to  derive  the  locations  from  the  arrival 
time  differences.  For  example,  in  the  method  of  Got  et  al.  (1994),  the  events  in  a  cluster  have 
precisely  the  same  take-off  angle  and  azimuth  to  each  station.  As  a  result,  the  derived  locations 
are  ultimately  relative,  not  absolute,  so  that  some  assumption  must  be  made  to  end  up  with 
useful  event  coordinates  (e.g.,  final  locations  are  computed  relative  to  a  catalog-based  cluster 
centroid).  Waldhauser  and  Ellsworth  (2000)  proposed  a  different  location  algorithm,  in  which  the 
spatial  partial  derivatives  for  the  set  of  events  are  evaluated  at  different  reference  points.  Thus 
the  absolute  event  locations  are  obtained,  rather  than  the  relative  locations.  It  is  assumed, 
however,  that  the  path  anomalies  from  velocity  heterogeneity  are  location  independent.  This 
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assumption  is  valid  for  closely  spaced  events,  but  is  not  true  for  far  apart  events.  As  a  result,  the 
event  locations  are  biased  due  to  velocity  heterogeneity  (Wolfe,  2003).  In  contrast,  the  latter 
approach  uses  the  relative  arrival  times  to  determine  a  much  smaller  number  of  adjusted  arrival 
time  picks,  but  these  picks  are  absolute  arrival  times  and  so  can  be  used  to  determine  absolute 
locations  (in  an  existing  velocity  model). 

We  have  developed  a  new  algorithm  that  combines  the  advantages  and  avoid  the 
disadvantages  of  the  above  approaches.  It  is  based  on  the  code  hypoDD  of  Waldhauser  and 
Ellsworth  (2000),  and  makes  use  of  both  absolute  and  relative  arrival  time  data.  The  method 
determines  a  three-dimensional  (3D)  velocity  model  jointly  with  the  absolute  and  relative  event 
locations.  This  approach  has  the  advantage  of  including  relative  arrival  times  with  their  quality 
values  along  with  absolute  arrival  times,  thereby  not  discarding  valuable  information  by  only 
using  adjusted  picks,  and  at  the  same  time  dispensing  with  simplifying  assumptions  about  ray 
path  geometries  or  path  anomalies  and  producing  absolute  locations,  not  just  relative  locations. 
Two  versions  have  been  developed,  one  Cartesian  and  the  other  spherical-earth. 

Methods,  Assumptions,  and  Procedures 

High-precision  waveform  cross-correlation 

WCC  can  be  a  powerful  tool  for  improving  earthquake  phase  identification  and 
subsequent  hypocenter  locations.  Our  software  is  being  built  upon  an  initial  package  that  was 
designed  to  retroactively  process  large  catalogues  having  existing  phase  picks  (P  and  S).  This 
package  compares  similar  waveforms  and  identifies  families  of  events  meeting  specified  cross¬ 
correlation  thresholds,  based  on  a  master  station  or  suite  of  stations.  Within  families  of  similar 
earthquakes,  waveforms  are  compared  on  a  station-by-station  basis  and  their  picks  are  adjusted 
by  applying  an  LI -norm  based  conjugate  gradient  technique  to  the  matrix  of  first  differences. 
Consistent,  mean-adjusted  picks  are  re-entered  into  trace  headers  for  other  uses  such  as 
waveform  stacking,  hypocenter  location  or  seismic  tomography. 

The  process  of  time  delay  calculation  and  bispectrum  verification  can  be  briefly  described 
as  follows.  For  the  k'th  mutual  station  that  provides  waveforms  for  two  earthquakes  i  and  j,  we 
rely  on  the  catalog  phase  (P  or  S)  picks  to  form  the  data  windows  used  for  time  delay 
calculations.  Then  using  the  band-pass  filtered  waveforms,  we  calculate  both  the  CC  value  CCij'‘ 
and  the  relative  time  delay  in  lag  number  for  the  event  pair.  If  the  value  of  CCjj'‘  is  larger 
than  a  threshold  CC*“'’,  we  also  perform  sub-sample  time  delay  calculation  to  get  by  a 

weighted  linear  fitting  of  the  cross  spectrum  phase  (Poupinet  et  al.,  1984).  After  we  make  the  CC 
calculations  for  all  the  mutual  stations  between  events  i  and  j,  we  can  obtain  the  maximum  CC 
value  CCij™*^  for  the  event  pair.  If  CCjj"’“’‘  is  larger  than  a  second  threshold  CC^®,  we  will  perform 
bispectrum  time  delay  estimation  for  each  mutual  station  of  the  event  pair.  We  obtain  two 
estimates  of  time  delay  in  lag  number,  i.e.  using  the  band-pass  filtered  waveforms  and 

A..k(bs2)  waveforms.  Thus  depending  on  the  extent  of  waveform  similarity  for  an 

event  pair,  which  is  measured  by  the  size  of  CCy™’',  four  possible  time  delay  estimations  are 
carried  out  at  a  mutual  station  for  either  P  or  S  waves. 
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Besides  the  aforementioned  threshold  value  that  researchers  often  use  to  select  time  delay 
estimates,  which  we  term  CC'""',  three  other  parameters  CC'""^  and  A*""  control  the 

bispectrum  verification  process.  Generally  a  CC'""'  a  CC'""^.  After  we  make  the  CC 
calculations  across  all  the  mutual  stations  of  an  event  pair,  we  compare  the  maximum  CC  value 
CCij"’®’‘  with  CC''"^.  If  CCij'"“  is  larger  than  CC*'"'^,  we  will  verify  the  CC  time  delay  estimates  for 
those  stations  with  CCy'^  larger  than  CC'""^.  A  CC  time  delay  estimate  is  trusted  (or  passes 
the  verification)  if  its  differences  from  both  Ay'‘^*“'^and  Ay'‘‘‘”^^,  the  ones  determined  with  the 
bispectrum  method,  do  not  exceed  the  tolerance  limit  A'"”.  In  other  words,  for  an  event  pair  with 
high  CCy'"“  value  we  will  select  the  CC  time  delay  estimates  as  long  as  they  pass  the  bispectrum 
verification,  even  though  the  CC  values  associated  with  some  of  them  are  smaller  than  CC'""*  and 
would  not  be  used  under  the  threshold  criterion  often  adopted  by  other  researchers.  If  CC''"’’  s 
Ccij™®*  <  cc'""^,  we  will  make  the  verifications  only  for  those  stations  with  CCy'‘  larger  than 
CC'""'.  If  CCy™’'  falls  below  CC''"’',  the  CC  time  delays  for  the  event  pair  are  simply  discarded. 
As  an  example,  we  apply  this  technique  to  obtain  bispectrum-verified  WCC  differential  times  for 
a  set  of  New  Zealand  earthquakes,  and  relocate  521  of  the  events  using  the  DD  algorithm 
hypoDD  (Waldhauser  and  Ellsworth,  2000).  We  find  that  the  bispectrum-verified  time  delays 
provide  more  clustered  earthquake  relocation  results  with  lower  arrival  time  residuals  compared 
to  the  threshold  criterion. 

We  do  not  claim  that  the  BS  method  always  obtains  a  better  time  delay  estimate  than  the 
CC  technique.  Rather,  we  use  the  BS  time  delay  estimates  to  check  the  reliability  of  the  ones 
computed  with  the  CC  technique.  By  applying  two  different  methods  that  work  in  different 
spectral  domains  (second-  and  third-order),  we  can  improve  the  reliability  of  the  selected  time 
delays. 

Generally  for  an  earthquake,  fewer  S  arrivals  are  available  in  the  phase  catalog  than  P 
arrivals  because  they  are  more  difficult  to  select.  We  can  obtain  additional  bispectrum-verified  S 
differential  times  for  those  waveforms  lacking  catalog  S  picks,  using  the  predicted  S  arrivals  either 
from  a  velocity  model  (Shearer,  1997)  or  from  the  S-P  time  of  a  nearby  event.  Our  work 
demonstrates  that  relocation  results  can  be  further  improved  with  these  additional  time 
constraints. 

Event  clustering 

We  tested  the  clustering  strategies  applied  with  the  agglomerative,  dendrogram-based 
clustering  algorithm  of  Rowe  et  al.  (2002a),  applying  variable  correlation  window  lengths  and 
targeting  different  parts  of  the  waveform  in  tiered  clustering  and  sub-clustering  processes.  We 
have  found  substantial  success  in  correlating  on  the  later  parts  of  a  waveform  rather  than  on  the 
first  arrival  in  some  instances,  as  signal-to-noise  of  the  first-P  arrival  can  be  low  for  emergent  or 
nodal  arrivals. 

Our  NMT  and  Sandia  collaborators  finalized  the  New  Mexico  GT  dataset  for  our  detailed 
clustering  and  real-time  tests.  The  New  Mexico  dataset  consists  of  two  parts,  one  spanning  July 
1 997-Februaty  1998  (Dataset  A)  that  has  undergone  a  complete  repicking,  relocation,  WCC,  and 


7 


clustering  analysis,  and  the  other  spanning  January  2001 -December  2002  (Dataset  B)  that  was 
prepared  for  real-time  testing. 

We  have  implemented  and  tested  a  two-step  clustering/location  technique  that  offers 
several  advantages  over  conventional  techniques  due  to  its  use  of  WCC  and  DD  location 
methods.  Step  1  event  clustering  is  based  solely  on  high-quality  GT  event  waveforms  without 
use  of  catalog  picks.  Step  2  location  estimates  offer  the  same  improvement  in  location  accuracy 
as  documented  by  the  use  of  high-resolution  DD  techniques  over  standard  techniques 
(Waldhauser  and  Ellsworth,  2000).  We  evaluated  the  technique  using  640  events  from  Northern 
California  (magnitude  range  4.3-5.3)  using  a  6°  x  6°  study  region  and  7  master  stations  from  the 
UC  Berkeley  seismic  network  (Figure  3).  In  step  1,  event  cluster  locations  are  estimated  via  CC 
of  Hilbert-envelope  waveforms.  A  test  event  is  cross-correlated  against  a  suite  of  reference  events 
binned  by  latitude,  longitude,  and  depth.  The  step  1  location  derived  from  the  CC  global 
maximum  is  then  used  as  a  starting  location  in  the  step  2  refinement  using  the  DD  relocation 
method  (hypoDD). 

We  initially  carried  out  testing  of  our  hierarchical  (multiple-frequency-band) 
clustering/WCC  analysis  on  a  California  GT  dataset  that  we  constructed  specifically  for  such 
testing  purposes.  The  dataset  consists  of  waveforms  for  697  events  with  magnitude  4  and  above 
from  central  and  northern  California.  We  affirmed  the  quality  of  the  GT  location  information  by 
carrying  out  a  DD  location  analysis.  Relocation  with  hypoDD  using  NCSN  data  from  25 
Berkeley  network  stations  provided  GT  hypocenter  information  for  the  test  events.  We  explored 
4  different  parameter  settings  controlling  the  extent  of  clustering  and  residual  outlier  rejection.  We 
found  that  the  DD  locations  reproduced  those  in  the  catalog  within  2  to  3  km  (in  epicenter  and 
depth)  on  average.  Then  the  dataset  was  divided  into  two  parts,  consisting  of  1 62  "master" 
events  (the  largest  in  each  location  bin)  and  the  remainder  as  "test"  events. 

An  adaptive,  automatic  phase-picking  and  epicenter-locating  program  is  now  fully 
developed  and  is  being  applied  to  update  the  New  Mexico  microseismicity  catalogue.  The 
programoperates  in  the  MATLAB  environment  (on  the  UNIX  platform)  under  MatSeis 
(developed  by  Sandia  National  Laboratories).  The  package  is  named  PLRR  (Pick-Locate-Repick- 
Relocate).  The  approach  of  PLRR  is  to  find  the  most  similar  event  from  a  reference  database,  and 
assign  phases  and  location  information  of  the  most  similar  event  in  a  reviewed  database  according 
to  the  spectrogram  cross-correlation.  No  preliminary  pick  or  location  information  is  needed  for 
the  unknown  event.  Some  characteristics  of  the  PLRR  package,  which  will  be  submitted  to  our 
Product  Integrator,  are: 

1 .  Data  format.  The  seismic  data  format  recognized  by  PLRR  is  the  CSS-3.0  flatfile  format. 

2.  Reference  database.  The  quality  of  phase-picking  and  epicenter-location  results  will 
ultimately  depend  on  the  coverage  and  quality  of  the  reference  database. 

3.  Spectrogram  waveform  cross-correlation.  We  find  the  phase  character  for  this  data  set  is  most 
diagnostic  and  discriminates  events  best  within  a  limited  bandwidth  between  6  and  1 0  Hz. 
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4.  Waveform-pair  match  and  event-pair  match.  We  have  explored  two  options  in  PLRR  to  find 
phases  and  locations  for  unknown  events.  The  first  one,  waveform-pair  matching,  finds  the  most 
similar  waveforms  in  the  reference  database  relative  to  all  waveforms,  station  by  station,  and 
assigns  the  corresponding  respective  phase  information  to  the  new  waveforms.  This  technique 
was  found  to  produce  too  many  false  associations.  The  second  one,  event-pair  matching,  is  more 
robust.  This  option  stacks  the  CC  coefficient  curves  of  waveform-pairs  between  the  unknown 
event  and  waveforms  grouped  from  each  event  in  the  reference  database,  and  finds  the  reference 
event  which  has  the  highest  stacked  CC  coefficient  value  (by  taking  the  mean/median  of  stacking 
CC  coefficient  curves)  as  the  most  similar  event  to  the  unknown  event. 

5.  Database  Operation.  Several  simple  functions  to  enhance  and  simplify  the  operation  of  CSS- 
3.0  database  are  added  into  PLRR,  such  as  the  "Update  monthly  table,"  which  add  new  event  to  a 
CSS-3.0  database. 

Waveform/wavelet  analyses 

We  developed  a  wavelet-based  auto-picker  for  detecting  and  picking  first-P  arrivals  that 
will  be  used  in  conjunction  with  our  location  processing  procedures,  both  for  identifying 
waveform  segments  for  correlation  and  for  adjusting  stack  picks.  A  manuscript  describing  the 
approach  is  in  the  final  stages  of  preparation.  Once  we  completed  the  construction  of  the 
California  GT  dataset,  we  began  to  test  potential  waveform  reshaping  methods  using  this  GT 
dataset.  Two  approaches  are  being  examined:  an  empirical  deconvolution  method  and  a 
deterministic  depth-based  method.  The  key  to  both  methods  is  to  determine  an  approach  that 
"simplifies"  the  first-arrival  wavelet  but  preserves  the  relative  arrival  time. 

We  tested  a  variety  of  decomposition  and  filtering  techniques  with  the  multiple  goals  of 
noise  reduction,  waveform  "homogenization"  (reshaping  to  similar  "wavelet"),  and  depth  phase 
identification.  Methods  examined  include  eigenimage  decomposition,  Karhunen-Loeve 
decomposition  (KLD),  Correlated  data  Adaptive  Noise  Canceling  (CANC),  and  coherency 
filtering.  A  key  issue  is  the  effect  of  arrival  misalignment  between  a  pair  of  traces  on  the  success 
of  a  given  technique.  We  tested  a  two-stage  process  whereby  the  similarity  and  strength  of  a 
CANC-filtered  trace  (relative  to  the  original  and  error  traces,  respectively)  is  used  to  align  the 
trace  to  the  reference  signal. 

Double-difference  tomography 

We  have  developed  a  new  DD  tomography  method  (Zhang  and  Thurber,  2003),  which 
takes  the  DD  location  method  (Waldhauser  and  Ellsworth,  2000)  a  step  further  to  solve  for  3D 
velocity  structure  simultaneously  with  hypocenter  locations  using  both  catalog  picks  and 
differential  arrival  times  (WCC  and  catalog).  At  the  local  scale  (10s  to  100s  of  kilometers),  the 
earth  can  be  represented  with  a  flat  model.  We  use  the  pseudo-bending  ray-tracing  algorithm  (Um 
and  Thurber,  1987)  to  find  the  rays  and  calculate  the  travel  times  between  events  and  stations. 
The  model  is  represented  as  a  regular  set  of  3D  nodes  and  the  velocity  values  are  interpolated  by 
using  the  tri-linear  interpolation  method.  The  hypocenter  partial  derivatives  are  calculated  from 
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the  direction  of  the  ray  and  the  local  velocity  at  the  source.  The  ray  path  is  divided  into  a  set  of 
segments  and  the  model  partial  derivatives  (calculated  in  terms  of  fractional  slowness 
perturbation,  so  that  the  derivatives  are  related  to  path  length)  are  evaluated  by  apportioning  the 
derivative  to  its  8  surrounding  nodes  according  to  their  interpolation  weights  on  the  segment 
midpoint  (Thurber,  1983). 

The  body  wave  arrival  time  T  from  an  earthquake  i  to  a  seismic  station  k  is  expressed 
using  ray  theory  as  a  path  integral 

T‘=T‘+/.'‘uds  (1) 


where  x*  is  the  origin  time  of  event  i,  u  is  the  slowness  field  and  ds  is  an  element  of  path  length. 
The  source  coordinates  (xi,  X2,  X3),  origin  times,  ray-paths,  and  the  slowness  field  are  unknown 
model  parameters.  The  relationship  between  the  arrival  time  and  the  event  location  is  nonlinear, 
so  a  truncated  Taylor  series  expansion  is  generally  used  to  linearize  Equation  (1).  The  misfit 
between  the  observed  and  the  predicted  arrival  times  can  be  linearly  related  to  the  perturbations 
to  the  hypocenter  and  velocity  structure  parameters 

1-1 


If  we  subtract  from  Equation  (2)  an  equivalent  equation  for  event),  we  have 


^  — ^Ax !  +  At '  +  J.'^bu  ds  -  ^  — ^Ax]  -  At|  -  J.'^bu  ds 
1.1  dx|  ‘  1.1  dx]  J 


(3) 


Assuming  that  these  two  events  are  near  each  other  so  that  the  paths  from  these  two  events  to 
the  common  station  are  almost  identical  and  the  velocity  structure  is  known,  then  Equation  (3) 
can  be  simplified  as 


■'  -rJ  = 


^dT'  •  dT-* 

’^^Axl  +  Atl-Y^AxJ-AtJ 
ax|  '  I'  iflaxj 


-Z 

1-1 


(4) 


Use  of  Equation  (4)  is  known  as  the  "double-difference"  (DD)  earthquake  location  algorithm 
(Waldhauser  and  Ellsworth,  2000).  Wolfe  (2003)  has  carried  out  an  elegant  comparison  of  this 
method  to  those  of  Jordan  and  Sverdrup  (1981)  and  Got  et  al.  (1994),  pointing  out  the  general 
similarities  and  subtle  differences,  and  discussing  the  advantages  and  limitations  of  each  method. 
Her  most  important  conclusion  is  "that  when  the  path  anomalies  from  velocity  heterogeneity 
change  strongly  with  earthquake  position,  the  bias  effects  can  be  reduced  in  the  relative  locations 
between  closely  spaced  earthquakes,  but  the  effects  cannot  be  reduced  in  the  relative  locations 
between  earthquakes  spaced  far  apart."  Two  assumptions  in  the  Waldhauser  and  Ellsworth 
(2000)  method  are  key:  an  a  priori  layered  velocity  model  is  assumed,  and  they  add  a  constraint 
equation  forcing  the  mean  shift  of  all  earthquakes  to  be  (approximately)  zero. 
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We  can  generalize  the  3D  DD  location  method  to  jointly  determine  the  3D  velocity 
structure  and  the  (absolute)  event  locations.  The  equations  for  the  DD  tomography  algorithm, 
tomoDD,  are: 


M  ‘  1=1  ax] 


(5) 


Our  tomography  algorithm  incorporates  first-difference  smoothing  constraints,  with  a  greater 
smoothing  constraint  weight  applied  to  horizontal  direction  than  the  vertical  (since  vertical 
velocity  changes  are  generally  expected  to  be  greater). 

At  the  regional  scale  (lOO's  to  lOOO's  of  km),  the  earth  should  not  be  treated  as  flat,  but  as 
a  sphere.  Some  velocity  discontinuities  such  as  Conrad,  Moho,  and  subducting  slab  boundary 
should  also  be  taken  into  account.  Finite-difference  ray  tracing  algorithms  developed  by  both 
Podvin  and  Lecomte  (1991)  and  Hole  and  Zelt  (1995)  are  able  to  accurately  calculate  first-arrival 
times  in  the  presence  of  extremely  severe,  arbitrarily  shaped  velocity  discontinuities.  These 
algorithms  solve  a  finite  difference  approximation  to  the  Eikonal  equation  in  the  regularly  gridded 
velocity  structure  through  a  systematic  application  of  Huygens'  principle.  This  procedure  can 
explicitly  take  into  account  different  propagation  modes  including  transmitted  and  diffracted 
body  waves,  and  head  waves  in  addition  to  direct  waves.  The  principle  of  reciprocity  is  utilized 
by  placing  the  source  point  at  the  seismic  station  location  and  interpolating  the  travel  times 
computed  at  the  grid  nodes  to  match  the  exact  earthquake  location  (Flanagan  et  al.,  1999). 

To  adapt  the  finite-difference  algorithm  to  the  regional  scale,  the  curvature  of  the  Earth  must  be 
taken  into  account.  Earth  flattening  transformation  is  only  valid  for  1 D  velocity  models  and  not 
for  3D  velocity  models.  Following  Flanagan  et  al.  (2000),  we  solved  this  problem  by 
parameterizing  a  spherical  surface  inside  the  Cartesian  volume  of  grid  nodes.  The  coordinate 
center  was  put  at  the  surface  of  the  Earth,  positive  X  and  Y-axes  point  to  the  direction  of  North 
and  West,  and  positive  Z-axis  points  downward,  respectively.  The  grid  nodes  above  the  Earth 
surface  (air  nodes)  are  attributed  to  the  true  velocity  values  traveling  in  the  air.  As  a  result,  all  the 
rays  travel  inside  the  Earth.  The  regular  uniform  velocity  (slowness)  grid  nodes  used  to  calculate 
the  travel  time  field  using  finite-difference  algorithms  are  interpolated  from  the  non-uniform 
inversion  grid  nodes  through  tri-linear  interpolation.  If  an  inversion  grid  node  is  2  km  above  the 
Earth's  surface,  then  it  is  treated  as  an  air  node  and  its  value  is  fixed  during  the  inversion.  We  treat 
each  station  as  a  source  and  calculate  travel  times  to  all  velocity  nodes  in  the  volume.  The  travel 
time  from  a  station  to  each  earthquake  is  interpolated  from  its  8  neighboring  nodes  through  tri- 
linear  interpolation.  The  ray  path  from  the  earthquake  to  the  station  is  found  iteratively  with 
increments  opposite  to  the  travel  time  gradient.  After  these  steps,  DD  tomography  is  adapted  to 
the  regional  scale  with  ray  paths  and  related  partial  derivatives  calculated  by  the  finite-difference 
method. 
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Results  and  Discussion 


High-precision  waveform  cross-correlation 

The  WCC  program  has  some  features  that  help  it  to  perform  in  an  auto-adaptive  manner. 
These  include:  adaptive,  full-spectrum  cross-coherency-based  prefiltering;  a  multi-component 
signal  covariance  estimator  that  provides  automatic,  adaptive  polarization  filtering  to  optimize 
signal  linearity;  and  the  use  of  multiple  eigentapers  in  the  subsample  cross-spectral  correlation 
step  to  reduce  signal  leakage  and  variance  when  estimating  the  signal  cross-spectra.  We  have 
adapted  the  cross-correlation  program  with  the  addition  of  eigentapers  in  the  coarse  (integer- 
sample)  cross-correlation  step,  and  we  are  in  the  process  of  comparing  this  approach  to  coarse 
lag  estimation  with  the  previous  method  of  applying  multiple  narrow-band  filtering.  Both 
approaches  provide  a  means  of  estimating  linearly  independent  realizations  of  the  time  series  and 
therefore  permitting  quantitative  error  estimates  for  the  integer  correlation  step;  in  future  work 
we  will  combine  the  two  strategies  via  application  of  multiple  eigenwavelet  methods. 

The  correlator  has  also  been  adapted  so  that  its  output  may  optionally  be  generated  in  a 
format  that  is  compatible  with  the  hypoDD  program,  so  that  joint  location  may  be  performed 
using  the  cross-correlation  comparisons  without  actually  solving  for  consistent  pick  lags  and 
correcting  the  picks  themselves.  Further,  the  correlator  now  has  an  option  of  generating  only  one 
set  of  first-differences,  for  applications  in  which  it  is  only  desired  to  compare  a  set  of  "reference 
events"  to  the  trace  of  a  single  ("test")  event.  This  option  is  being  tested  in  our  hierarchical 
approach  described  in  the  next  section. 

Figure  1  demonstrates  the  application  of  the  CC  and  bispectrum  methods  to  the 
waveforms  of  two  closely  spaced  M2.5  New  Zealand  earthquakes.  The  waveforms  for  these 
events  recorded  at  nearby  stations  are  highly  similar,  but  the  signal  to  noise  ratios  (SNRs)  for  the 
records  at  station  MOW  (epicentral  distance  of  93  km)  are  noticeably  lower  (Figs.  1 A  and  1 B). 
The  two  bispectrum-correlation  series  (Figs.  IE  and  IF)  peak  at  lag  6  and  7  respectively,  whereas 
the  two  CC  series  (Figs.  1C  and  ID)  reach  maxima  at  lags  very  different  from  each  other  and  also 
from  the  two  bispectrum-correlation  series.  Careful  examination  reveals  that  the  two  CC  series 
reach  local  maxima  at  lag  7.  They  fail  to  become  global  maxima  there,  however,  apparently  due  to 
the  contamination  of  correlated  noise.  Figures  IG  and  IH  plot  the  stacked  signals  for  both  raw 
and  band-pass  filtered  waveforms.  In  both  cases  the  lower  trace  (the  stacked  signal  obtained  by 
shifting  the  waveform  of  the  first  event  with  bispectrum-determined  lags)  has  a  larger  root-mean- 
square  (RMS)  amplitude  than  the  upper  trace,  which  is  obtained  after  shifting  the  waveform  of 
the  first  event  with  the  CC-determined  lags.  The  above  example  corroborates  the  findings  of 
Nikias  and  Pan  (1988)  and  Yung  and  Ikelle  (1987)  that  the  bispectrum  method  works  better  than 
CC  techniques  in  getting  relative  time  delay  between  two  similar  signals  contaminated  by 
correlated  noise. 

Figure  2  shows  examples  of  the  application  of  the  BS  verification  method  to  waveforms 
recorded  at  New  Zealand  stations  for  two  pairs  of  relatively  close  events.  In  Figure  2a,  the  CC 
value  (0.74)  is  above  the  nominal  threshold  (0.70),  but  the  lag  does  not  agree  with  that  obtained 
from  BS  (bottom  panel),  so  this  CC  value  is  rejected.  From  the  seismograms,  we  can  see  that  this 
is  apparently  a  case  of  cycle  skipping  -  the  filtered  waveforms  coincidentally  correlate  slightly 
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better  at  a  lag  of  -10.  In  Figure  2b,  both  methods  provide  the  same  time  delay,  but  the  associated 
CC  coefficient  is  only  0.50.  Under  the  threshold  criterion,  such  a  time  delay  estimate  associated 
with  a  low  CC  coefficient  is  simply  discarded.  Using  the  BS  checking  method,  however,  we  can 
identify  the  reliable  time  delays  from  these  "seemingly  not  very  similar"  waveforms,  such  as  for 
this  pair,  and  as  a  result  provide  more  control  over  the  relative  locations  of  the  events.  Du  et  al. 
(2004)  apply  this  technique  to  822  New  Zealand  earthquakes  in  the  Wellington  area  and  find  that 
the  CC  time  delays  verified  with  the  BS  method  provide  improved  earthquake  relocation  results 
(noticeably  more  clustered)  compared  to  those  selected  with  the  standard  threshold  criterion.  An 
example  of  the  improvement  of  the  BS-verified  locations  compared  to  original  catalog  locations, 
DD  locations  using  differential  catalog  picks,  and  DD  locations  using  standard  threshold  CC  lags 
is  shown  in  Figure  3  for  a  sub-region  near  Lake  Wairarapa,  NZ. 

Event  clustering 

Waveform  clustering,  used  to  segregate  signals  in  to  families  of  similar  events,  has  been 
tested  with  a  variety  of  approaches.  One  successful  approach  has  been  to  identify  clusters  based 
on  very  high  standards  of  cross-correlation  (0.8  or  0.85),  resulting  in  many  small  multiplets  of 
highly  similar  events.  These  are  then  cross-correlated  and  consistent  pick  adjustments  are 
determined  via  a  conjugate  gradient  solver.  Waveforms  are  aligned  on  their  new  picks  and  stacked 
to  provide  a  composite  waveform  representation  for  the  cluster.  These  stacks  are  then  cross- 
correlated  among  the  multiplets  and  inconsistencies  in  the  mean  adjusted  picks  can  be  corrected 
by  finding  consistent  adjustments  among  the  stacks.  The  resulting  stack  pick  lags  are  propagated 
to  the  parent  traces.  This  has  been  shown  to  reduce  the  inter-cluster  inconsistencies  that  result  in 
cluster  centroid  biases  when  earthquakes  are  located  using  the  new  picks  (these  biases  existed  in 
the  initial  catalogue  as  well).  Another  related  approach  that  has  shown  good  success  is 
application  of  the  clustering  program  to  the  waveform  stacks.  Those  traces  whose  stacks  cluster 
together  based  on  a  high  similarity  coefficient  are  joined  into  larger,  second  order  families.  Their 
member  traces  are  cross-correlated  and  adjusted  for  consistency  within  this  larger  group.  The 
second-order  clusters,  following  pick  adjustment,  are  stacked  yet  again  and  their  stacks  are 
compared  for  similarity.  This  process  may  be  performed  in  a  hierarchical  manner  until  stacks  fail 
to  cluster.  Final  adjustment  of  mean  centroid  picks  on  the  ultimate  stacks  has  been  tested  using 
an  autopicker;  this  works  well  when  many  traces  contribute  to  the  stack  but  we  are  still  testing 
different  methods  of  estimating  pick  uncertainty  from  the  autopicker  functions.  Phase-weighted 
stacking,  jackknifing  methods  or  multiple-frequency-band  approaches  may  provide  the  best 
uncertainty  estimates. 

We  constructed  two  regional  ground-truth  datasets,  one  for  New  Mexico  and  one  for 
California.  The  former  dataset  is  from  the  New  Mexico  Tech  short-period  network  and  is 
dominated  by  small  explosions,  whereas  the  latter  is  from  the  UC-Berkeley  broadband  network 
and  contains  relatively  large  earthquakes.  The  New  Mexico  dataset  consists  of  two  parts,  one 
spanning  July  1 997-February  1998  that  has  undergone  a  complete  repicking,  relocation,  WCC, 
and  clustering  analysis,  and  the  other  spanning  January  2001 -present  that  is  being  prepared  for 
real-time  testing.  The  California  dataset  was  extracted  from  the  Northern  California  Earthquake 
Data  Center  (NCEDC).  We  first  affirmed  the  quality  of  the  NCEDC  catalog  location  information 
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by  carrying  out  a  DD  location  analysis  using  catalog  time-difference  data  for  1 30  events  and  23 
stations.  We  found  that  the  DD  locations  reproduced  those  in  the  catalog  within  2  to  3  km  (in 
epicenter  and  depth)  on  average.  We  then  used  a  combination  of  automatic  and  manual  picking  to 
augment  and  improve  the  available  P  arrival  picks,  and  to  provide  analyst  pick  weights. 

We  present  the  results  of  hierarchical  WCC/clustering  analysis  on  the  California  GT 
dataset,  based  on  WCC  analysis  of  P  through  surface  waves  in  a  first  step  and  the  P  phase  in  a 
second  step.  We  binned  the  events  in  0.1°  latitude  and  longitude  bins  and  in  2.5  km  depth  bins. 
The  largest  magnitude  events  per  bin  comprise  our  1 62  reference  events,  and  the  next  largest  24 
events  are  used  for  the  test  dataset.  The  reference  events  are  displayed  in  Figure  4a,  and  the  test 
events  and  stations  used  for  the  preliminary  hierarchical  analysis  are  shown  in  Figure  4b.  The 
WCC  analysis  for  the  two  steps  was  accomplished  in  a  similar  manner.  Waveform  data  from  each 
test  event  were  cross-correlated  against  the  reference  event  waveforms  to  determine  the  maximum 
value,  in  the  whole  region  using  low-frequency  data  in  step  one,  and  in  a  2°  region  about  the  step- 
one  location  estimate  using  high-frequency  data  in  step  two.  The  maximum  cross-correlation 
value  per  bin  from  the  7  stations  was  contoured  as  an  indicator  of  test  event  location.  Figures  5a 
and  b  display  the  results  for  a  test  event  located  near  Los  Angeles.  For  this  example,  the  WCC 
analysis  is  performed  on  the  vertical  component  data  for  both  steps.  The  estimated  location  for 
this  event  falls  with  25  km  of  the  catalog  location  using  the  low-frequency  data,  and  improves  to 
7.5  km  using  high-frequency  data. 

We  then  tested  the  two-step  event  location  technique  with  a  California  GT  data  set.  We 
augmented  the  existing  P  and  S  picks  from  the  network  catalog  using  skew  and  kurtosis  analysis 
followed  by  a  manual  evaluation.  Many  catalog  P  picks  required  manual  adjustment.  We  then 
prepared  WCC  data  based  on  3-  and  5-second  P  windows  (Rowe  et  al.,  2002a).  In  step  1,  event 
locations  are  estimated  via  WCC  of  Hilbert-envelope  waveforms.  The  entire  seismogram  was 
filtered  and  an  eigenvector  decomposition  of  the  covariance  matrix  was  applied  to  the  3- 
component  data  (Rowe  et  al.,  2002a)  prior  to  the  calculation  of  Hilbert-envelopes  for  the 
principal  direction.  Test  event  data  for  each  of  the  master  stations  were  cross-correlated  against 
1 6 1  GT  events  with  the  position  of  the  CC  global  maximum  determining  the  step  1  location 
(Figure  3).  The  Step  1  procedure  accurately  located  the  epicenters  within  20  km  (approximately 
1  bin)  of  the  NCEDC  catalog  location  for  1 6  of  the  20  test  events.  However,  depth  estimates 
were  poorly  constrained.  The  Step  2  refinement  based  on  WCC  of  the  high-frequency  P  arrivals 
followed  by  DD  relocation  improved  both  epicenter  and  depth  estimates.  For  each  test  event, 
subsets  of  the  catalog  and  WCC  data  are  prepared  for  use  in  hypoDD  with  the  step  1  location 
providing  the  initial  test-event  location.  This  two-step  procedure  offers  advantages  over 
conventional  location  techniques  including  a  Step  1  event  location  not  dependent  on  the  accuracy 
of  catalog  picks  and  a  Step  2  improvement  in  location  accuracy  using  the  high-accuracy  DD 
technique.  We  submitted  the  automated  two-step  location  code  (PLRR  UW)  to  our  Product 
Integrator. 

Waveform/wavelet  analyses 

We  tested  a  variety  of  decomposition  and  filtering  techniques  with  the  multiple  goals  of 
noise  reduction,  waveform  "homogenization"  (reshaping  to  similar  "wavelet"),  and  depth  phase 
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identification.  Methods  examined  include  eigenimage  decomposition,  Karhunen-Loeve 
decomposition  (KLD),  Correlated  data  Adaptive  Noise  Canceling  (CANC),  and  coherency 
filtering.  A  key  issue  is  the  effect  of  arrival  misalignment  between  a  pair  of  traces  on  the  success 
of  a  given  technique.  Some  methods  rely  on  signals  being  time-aligned  before  processing,  which  of 
course  cannot  be  assumed  for  seismic  data  in  a  monitoring  context.  We  also  developed  a  wavelet- 
based  auto-picker  for  detecting  and  picking  first-P  arrivals  that  can  be  used  in  conjunction  with 
our  location  processing  procedures,  both  for  identifying  waveform  segments  for  correlation  and 
for  adjusting  stack  picks  (Zhang  et  al.,  2003). 

We  have  been  testing  the  utility  of  the  KLD  method  for  extracting  similar  signals  from 
noisy  data.  A  synthetic  example  is  shown  in  Figure  6.  In  this  particular  case,  we  construct 
initially  noise-free  seismograms  (upper  left  panel)  that  consist  of  a  first-arriving  signal  that  is 
time-aligned  and  an  identical  second-arriving  signal  with  substantial  moveout.  Gaussian  noise 
with  RMS  amplitude  20%  of  the  signal  amplitude  is  added  (upper  right  panel)  and  eigenimage 
decomposition  (not  shown)  and  KLD  are  carried.  The  eigenimage  decomposition,  using 
eigenvalues  of  average  value  and  greater,  preserves  the  time-aligned  signal  but  degrades  the 
secondary  arrival,  leaving  most  of  the  latter's  energy  in  the  residual  "noise."  KLD  (middle  left 
panel),  again  using  eigenvalues  of  average  value  and  greater  (bottom  left  panel),  preserves  both  the 
time-aligned  signal  and  the  secondary  arrival  (middle  right  panel). 

For  a  second  example,  we  apply  KLD  to  real  data  from  a  group  of  earthquakes  of 
magnitude  3.6  to  3.8  in  Long  Valley,  CA,  from  January  1998.  In  the  catalog,  the  events  are  within 
about  1  km  of  each  other  in  epicenter  but  have  estimated  depths  ranging  from  0.5  to  5.0  km. 
Relocation  using  hypoDD  and  catalog  arrival  times  plus  catalog-generated  arrival  time  differences 
yields  epicenters  within  less  than  1  km  and  depths  within  250  m  of  each  other.  The  analyzed 
waveforms  were  recorded  at  station  ORV  with  an  epicentral  distance  of  314  km. 

The  results  of  the  KLD  analysis  are  shown  in  Figure  7.  The  raw  data  traces  are  relatively 
noisy.  In  this  case,  we  recognize  that  the  eigenvectors  for  the  four  largest  eigenvalues  represent 
long-period  noise  (bottom  left  panel),  and  we  exclude  them  along  with  the  smaller  eigenvalues. 

The  resulting  KLD-processed  traces  (upper  right  panel)  are  quite  similar  and  the  data  residue 
(middle  left  panel)  contains  no  obvious  remnant  signal  of  interest.  In  contrast,  the  eigenimage 
decomposition  does  little  if  anything  to  improve  the  signal  similarity,  only  removing  some  of  the 
higher  frequency  noise. 

An  alternative  approach  to  noise  reduction  via  extracting  common  waveform  signals  is 
adaptive  filtering.  We  have  experimented  with  a  CANC  algorithm  of  Hattingh  (1988).  Input  into 
the  CANC  procedure  consists  of  a  "reference  signal"  SI  and  a  "primary  signal"  S2,  correlated  to 
each  other  in  some  unknown  way.  Each  signal  is  presumed  to  contain  noise  that  is  uncorrelated 
with  the  noise  from  the  other  trace.  The  reference  signal  is  used  as  a  "learning  signal"  to  derive 
filter  coefficients.  The  filter  adjusts  its  impulse  response  to  minimize  the  mean-squared  error 
between  filter  output  and  the  primary  signal  (in  terms  of  power).  The  algorithm  rapidly  solves 
for  the  filter  coefficients  without  requiring  matrix  inversion  or  a  CC  operation.  The  CANC 
procedure  is  an  iterative  solution,  computing  updated  filter  coefficients  at  each  iteration. 

We  apply  the  CANC  procedure  to  seismic  data  from  two  events  in  a  cluster  with  similar 
catalog  epicenters  but  different  depths,  recorded  at  UC-Berkeley  network  station  SAG.  The  two 
events  are  from  the  1994  Northridge  aftershock  sequence  and  are  located  within  7  km  of  each 
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other  in  epicenter  based  on  catalog  information.  Trace  1  (Figure  8a),  the  reference  signal,  is  from  a 
shallow  earthquake  (2  km  depth),  whereas  trace  2  (Figure  8b),  the  primary  signal,  is  from  an 
earthquake  occurring  at  greater  depth  (12  km).  We  apply  adaptive  filtering  to  shape  the 
windowed  signal  from  the  deeper  earthquake  to  extract  the  signal  common  to  the  shallow  event  in 
an  attempt  to  improve  subsequent  WCC  results.  We  use  a  7  s  window  containing  approximately 
4.5  s  of  data  following  the  P  pick.  The  windowed  signal  from  the  deeper  event  (Figure  8b)  was 
aligned  relative  to  the  shallow  event  signal  (Figure  8a)  initially  by  the  P  picks  and  then  fine-tuned 
by  applying  a  lag  search.  The  zero-lag  CC  between  the  primary  and  filtered  traces  (Figure  8c)  and 
the  signal-to-noise  ratio  of  the  filtered  trace  relative  to  the  error  trace  (Figure  8d)  have  coincident 
maxima.  After  5  iterations  of  the  CANC  algorithm,  the  signal  from  the  deeper  event  is 
successfully  decomposed  into  an  estimate  of  the  common  signal  (Figure  8e)  and  the  dissimilar 
error  signal  (Figure  8f). 

We  can  compare  the  estimate  of  common  signal  to  individual  trace  output  by  the 
application  of  coherency  filtering.  Shown  are  the  coherency-filtered  traces  for  the  shallow  and 
deep  event  (Figures  8g  and  h).  The  frequency  content  of  the  two  clustered  events  is  very  similar 
as  might  be  expected.  Consequently,  coherency  filtering  yields  output  traces  similar  to  the 
original  trace  data.  The  adaptive  filter  shaping  of  the  deep  event's  waveform  visually  yields  a 
signal  (Figure  8e)  closer  in  appearance  to  the  shallow  event  (Figure  8a).  Application  of  adaptive 
filtering  increases  the  maximum  CC  value  for  the  shallow  event  with  the  CANC-filtered  trace 
from  0.5 1  to  0.83.  CC  of  the  coherency-filtered  data  yields  a  maximum  value  of  0.70.  These  CC 
maxima  occur  at  the  same  lag.  Thus  the  CANC  approach  can  successfully  extract  a  more  similar 
waveform  from  the  primary  signal  than  coherency  filtering  while  preserving  timing. 

Double-difference  tomography 

We  have  tested  the  local-scale  version  of  DD  tomography  (tomoDD)  with  a  synthetic 
dataset  that  was  constructed  on  the  basis  of  an  idealized  model  of  the  velocity  structure  of  the 
San  Andreas  Fault  in  central  California  (Kissling  et  al.,  1994)  shown  in  Figure  9.  The  events  and 
stations  used  to  construct  the  synthetic  dataset  are  from  the  actual  seismicity  and  USGS  stations 
in  the  Loma  Prieta  region  (Figure  7).  We  added  Gaussian  random  noise  with  zero  mean  and  a 
standard  deviation  of  0.04  s  to  all  the  true  arrival  times.  In  addition,  we  also  added  a  constant 
noise  term  to  the  arrivals  at  each  station  from  a  uniform  distribution  between  -0.3  s  and  0.3  s. 
This  simulates  the  case  that  the  systematic  errors  (model  errors  and  pick  bias)  associated  with 
the  arrival  times  are  larger  than  the  random  ones.  We  construct  the  pseudo  "cross-correlation" 
data  directly  from  the  absolute  arrival  times  by  differencing  the  synthetic  arrival  times  at  common 
stations  for  events  pairs  within  20  km.  As  a  result,  the  "cross-correlation"  data  are  more  accurate 
than  the  absolute  data. 

We  first  use  the  DD  location  algorithm  hypoDD  to  relocate  the  events  (Waldhauser, 
2001).  The  ID  velocity  model  used  is  the  Dietz  and  Ellsworth  (1990)  model,  based  on  seismic 
refraction  and  earthquake  data.  We  use  both  the  residual  weighting  and  distance  weighting 
schemes  in  the  inversion  to  control  the  large  residuals,  as  did  Waldhauser  and  Ellsworth  (2000). 
The  absolute  differences  between  this  set  of  event  relocations  and  the  true  locations  and  their 
standard  deviations  are  shown  in  Table  1 .  The  results  indicate  that  the  event  relocations  from  the 
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DD  location  algorithm  based  on  a  ID  velocity  model  have  a  substantial  bias  (>1  km  in  each 
coordinate  direction)  from  the  true  locations.  This  bias  is  caused  by  a  difference  between  the 
velocity  model  (horizontal  layers)  used  to  calculate  the  DD  locations  and  the  true  velocity  model 
(vertical  "sandwich")  used  to  generate  the  data.  The  heterogeneity  of  the  true  velocity  model 
makes  path  anomalies  different  for  different  events.  However,  the  DD  location  method  assumes 
that  the  path  anomalies  are  location-independent,  and  this  assumption  introduces  bias  into  event 
locations  (Wolfe,  2003). 

Next,  we  show  the  results  of  applying  the  standard  tomography  method  to  the  relatively 
noisy  absolute  arrival  times.  The  inversion  starts  from  the  same  ID  velocity  model  as  the  DD 
location  method.  The  X-Y  nodes  used  to  represent  the  velocity  model  are  shown  in  Figure  10;  in 
depth,  nodes  are  positioned  at  0, 3,  7,  11,  and  16  km.  The  computational  algorithm  is  identical  to 
that  for  DD  tomography,  but  the  differential  times  are  excluded.  Both  damping  and  smoothing 
(with  weight  5.0  in  both  the  horizontal  and  vertical  directions)  are  applied  to  the  inversion  to 
make  the  solution  more  stable.  Figure  8a  shows  horizontal  slices  through  the  velocity  structure 
obtained  from  the  standard  tomography.  The  absolute  difference  between  this  velocity  model  and 
the  true  velocity  model  has  a  median  value  of  0. 1 64  km/s,  a  mean  value  of  0.245  km/s  and  a 
standard  deviation  of  0.249  km/s.  The  main  features  of  the  true  velocity  model  are  evident  in  the 
depth  slices.  The  event  locations  have  smaller  errors  than  those  from  the  DD  location  method 
(Table  1). 

DD  tomography  uses  both  the  noisy  absolute  arrival  times  and  more  accurate  differential 
times.  We  also  apply  damping  and  the  smoothing  constraints  to  stabilize  the  solution,  with  the 
residual  weighting  used  to  penalize  the  large  residuals  during  the  inversion.  DD  tomography 
starts  with  the  same  1 D  velocity  model  and  uses  the  same  inversion  grid  as  the  standard 
tomography.  Figure  1 1  a  and  b  shows  the  horizontal  slices  of  the  velocity  structure  from 
conventional  and  DD  tomography.  The  absolute  difference  between  the  DD  tomography  velocity 
structure  and  the  true  velocity  structure  has  a  median  value  of  0. 136  km/s,  a  mean  value  of  0. 1 78 
km/s,  and  a  standard  deviation  of  0. 1 64  km/s,  compared  to  a  median  value  of  0. 1 64  km/s,  a  mean 
value  of  0.245  km/s,  and  a  standard  deviation  of  0.249  km/s  for  conventional  tomography.  DD 
tomography  characterizes  well  the  low  velocity  zone  of  the  true  velocity  model,  especially  at  the 
depths  of  3  km  and  7  km.  Subtracting  the  DD  tomography  solution  and  the  standard  tomography 
solution  from  the  true  model  (Figure  1  Ic  and  d),  we  find  that  the  velocity  model  from  DD 
tomography  has  a  more  correct  value  in  the  low-velocity  trough  except  at  the  depth  of  0  km, 
where  the  ray  paths  from  event  pairs  almost  completely  overlap  and  the  accuracy  of  velocity 
structure  is  mainly  controlled  by  absolute  catalog  data.  This  indicates  that  DD  tomography 
recovers  better  the  low-velocity  zone,  and  overall  the  velocity  model  is  more  accurate  than  the 
standard  tomography.  Compared  with  the  standard  tomography  method,  DD  tomography  also 
produces  more  accurate  event  locations  (Table  1). 

Conclusions  and  Recommendations 

Advances  in  techniques  for  waveform  alignment,  event  clustering,  waveform 
decomposition,  and  location  and  tomography  have  been  achieved  with  the  support  of  this 
contract.  We  echo  the  conclusion  voiced  by  Prof  Paul  Richards  that  a  quantum  change  is  needed 
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in  the  manner  in  which  earthquake  locations  are  determined  in  the  seismic  monitoring 
environment  (both  nuclear  and  otherwise).  Effective  tools  for  event  clustering,  waveform 
alignment,  and  joint  location/tomography  are  available,  and  a  concerted  effort  will  be  required  to 
implement  these  tools  in  "real-world"  monitoring  situations.  Accomplishing  this  will  probably 
require  closer  and  more  direct  interaction  between  the  tool  developers  and  the  users,  in  particular 
to  allow  for  direct  feedback  from  the  users  to  the  developers. 

Differential  arrival  times  for  pairs  of  seismic  events  observed  at  the  same  station  are  often 
calculated  by  WCC  techniques.  Researchers  often  choose  the  differential  times  to  use  based  on 
the  associated  cross-correlation  (CC)  values  exceeding  a  specified  threshold.  When  two  similar 
time  series  are  corrupted  by  correlated  noise,  the  time  delay  estimate  calculated  with  the  WCC 
technique  may  not  be  reliable.  Thus  the  selection  of  a  threshold  value  is  important.  If  it  is  set  too 
high,  then  only  a  limited  number  of  very  accurate  differential  time  data  are  available  to  constrain 
the  relative  positions  of  earthquakes.  If  the  threshold  value  is  set  too  low,  then  many  unreliable 
differential  time  estimates  are  used  and  they  will  negatively  affect  the  final  relocation  results.  The 
bispectrum  method  can  suppress  the  correlated  Gaussian  noise  sources  in  two  similar  time  series 
and  be  used  to  obtain  the  relative  time  delay  between  them.  We  compute  time  delay  estimates 
between  the  waveform  pairs  with  the  bispectrum  method  and  use  them  to  verify  the  WCC- 
determined  one.  This  technique  can  provide  quality  control  over  the  selected  time  delay  estimates 
and  potentially  provide  more  differential  travel  time  data  for  close  events  pairs,  by  verifying  the 
reliability  of  differential  times  that  might  not  meet  the  threshold  criterion. 

WCC  is  vulnerable  to  problems  associated  with  correlated  noise  that  can  be  overcome 
using  bispectrum  analysis.  It  is  our  recommendation  that  a  verification  approach  comparable  to 
that  described  here  be  adopted  for  earthquake  location  applications  for  which  very  high  accuracy 
is  essential. 

A  hierarchical  location  approach,  for  "new"  events  relative  to  a  set  of  reference  events, 
that  starts  with  CC  of  long  records  at  low  frequencies  to  identify  the  approximate  source  region 
based  solely  on  waveform  or  spectrogram  similarity  and  then  switches  to  single-event  or  DD 
location  using  CC  of  short  records  at  high  frequencies  works  effectively  for  our  regional  datasets. 
We  have  tested  two  versions  of  a  two-step  event  location  technique  with  California  and  New 
Mexico  GT  data  set.  In  step  1  of  the  first  (UW)  version,  event  locations  are  estimated  via  WCC 
of  Hilbert-envelope  waveforms.  The  Step  1  procedure  accurately  located  the  epicenters  within  20 
km  (approximately  1  bin)  of  the  NCSN  catalog  location  for  16  of  the  20  test  events.  However, 
depth  estimates  were  poorly  constrained.  The  Step  2  refinement  based  on  WCC  of  the  high- 
frequency  P  arrivals  followed  by  DD  relocation  improved  both  epicenter  and  depth  estimates. 
This  two-step  procedure  offers  advantages  over  conventional  location  techniques  including  a  Step 
1  event  location  not  dependent  on  the  accuracy  of  catalog  picks  and  a  Step  2  improvement  in 
location  accuracy  using  the  high-accuracy  DD  technique.  For  the  other  (NMT)  version,  the 
approach  is  to  find  the  most  similar  event  from  a  reference  database,  and  assign  phases  and 
location  information  of  the  most  similar  event  in  a  reviewed  database  according  to  the 
spectrogram  cross-correlation.  No  preliminary  pick  or  location  information  is  needed  for  the 
unknown  event.  Both  approaches  have  proven  effective  in  our  trials.  This  approach  should  next 
be  tested  in  a  real-time  monitoring  environment. 
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We  tested  a  variety  of  decomposition  and  filtering  techniques  with  the  multiple  goals  of 
noise  reduction,  waveform  "homogenization"  (reshaping  to  similar  "wavelet"),  and  depth  phase 
identification.  The  KLD  and  CANC  methods  show  the  most  promise  for  allowing  the  isolation  of 
similar  underlying  waveforms  that  are  not  time-aligned. 

We  have  developed  the  new  method  of  DD  tomography.  The  inversion  code  tomoDD 
takes  the  double-difference  location  method  (hypoDD)  a  step  further  to  solve  for  3D  velocity 
structure  simultaneously  with  hypocenter  locations  using  both  catalog  picks  and  differential 
arrival  times  (WCC  and  catalog).  Compared  to  conventional  tomography,  the  result  is  a 
sharpening  of  the  seismicity  distribution  equal  in  quality  to  hypoDD  plus  a  sharpening  of  the 
velocity  image  due  to  removal  of  the  majority  of  the  location  scatter.  We  have  developed  local- 
and  regional-scale  versions  of  tomoDD,  and  have  applied  each  to  synthetic  and  real  datasets. 
Synthetic  tests  show  the  method  yields  locations  with  greater  accuracy  than  either  DD  location 
or  standard  tomography,  and  a  velocity  model  that  is  closer  to  the  true  model  than  standard 
tomography.  The  DD  tomography  approach  works  with  either  waveform-based  or  catalog-based 
differential  time  data,  and  thus  is  applicable  to  any  dataset  to  which  standard  tomography  can  be 
applied. 
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Figure  1.  (A)  Windowed  waveforms  (1  second  before  and  2.82  seconds  after  the  catalog  P  pick) 
for  two  New  Zealand  earthquakes  recorded  at  station  MOW.  The  station  sampling  rate  is  50 
sample/sec  and  the  signals  are  aligned  by  the  P  picks.  (B)  Windowed  waveforms  band-pass 
filtered  between  1  and  6  Hz.  (C,  D)  CC  series  for  the  (C)  unfiltered  (raw)  and  (D)  filtered 
waveforms.  (E,F)  Bispectrum-correlation  series  for  the  (E)  raw  and  (F)  filtered  waveforms.  (G) 
Stacked  raw  waveforms.  The  upper  trace  is  obtained  by  shifting  the  waveform  of  the  first  event 
with  the  CC  determined  lag  relative  to  the  second  one,  while  the  lower  trace  is  computed  after 
shifting  the  waveform  of  the  first  event  with  the  lag  calculated  with  the  bispectrum  method.  The 
root-mean-square  amplitudes  for  the  two  traces  are  0.25  and  0.27  respectively.  (H)  Same  as  G  for 
the  stacked  band-pass  filtered  waveforms.  The  root-mean-square  amplitudes  for  the  two  traces 
are  0.32  and  0.34  respectively. 
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Windowed  Waveforms  (KIW;  50  sp/sec;  filtered)  Windowed  Waveforms  (OTW;  50  sp/sec;  filtered) 


Figure  2.  Examples  of  a)  rejection  of  a  high  CC  value  due  to  inconsistent  results  between  the  CC 
and  BS  lags  and  b)  acceptance  of  a  low  CC  value  due  to  consistent  results  between  the  CC  and 
BS  lags. 
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Figure  3.  Locations  of  53  earthquakes  near  lake  Wairarapa,  NZ.  (A)  Before  relocation;  (B) 
relocated  with  catalog  differential  travel  times;  (C)  relocated  with  CC  differential  travel  times 
chosen  with  the  threshold  criterion;  (D)  relocated  with  CC  differential  travel  times  verified  with 
the  bispectrum  method;  (E)  relocated  with  both  catalog  and  threshold-chosen  CC  differential 
travel  times;  (F)  relocated  with  both  catalog  and  bispectrum-verified  CC  differential  travel  times. 
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Figure  4.  Maps  of  (a)  reference  events  and  (b)  test  events  and  test  stations  from  our  California 
ground-truth  dataset. 


Figure  5.  Contour  plots  (map  view  and  depth  sections)  of  the  maximum  cross-correlation  value 
per  bin  for  the  7  stations  for  a  test  event  (located  at  lat.  34.3°,  Ion.  -1 18.58°,  and  depth  1.8  km, 
cross)  correlated  against  the  set  of  reference  events  in  the  (a)  low-frequency  band  (step  one)  and 
(b)  high-frequency  band  (step  two).  Note  the  good  low-resolution  location  estimate  on  the  left 
(circle)  and  the  improved  high-resolution  location  estimate  (circle)  on  the  right. 
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Synthetic  Data  with  Moveout  of  1 .0  seconds  per  trace 


Raw  Data  with  gaussian  noise  (0.20)  added 
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Figure  6.  Example  application  of  KLD  to  a  set  of  seismograms  with  a  pair  of  "arrivals"  only  one 
of  which  is  time-aligned  across  the  records  (top  left),  with  noise  added  (top  right).  The  data 
reconstructed  from  the  61  eigenvectors  with  eigenvalues  above  average  in  size  (bottom  left) 
retains  both  the  time-aligned  arrival  and  the  arrival  with  moveout  (middle  left),  leaving  little  signal 
energy  in  the  residual  image  (middle  right).  In  contrast,  eigenimage  decomposition  cannot 
reconstruct  both  arrivals  adequately. 
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Raw  Data  (4  traces;  481  points  per  trace) 


Constructed  Data  with  #5  to  #34  eigenvectors 
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Data  Residue 


First  50  eigenvalues 


Plot  of  the  first  15  eigenvectors 
fvWWWWWWWWq 


Plot  of  reconstmction  coefficients 
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Figure  7.  Example  application  of  KLD  to  real  data  (upper  left).  Large  and  small  eigenvalues  were 
excluded  in  this  case.  The  reconstructed  data  (upper  right)  from  the  30  eigenvectors  with 
eigenvalues  from  #5  to  #34  (middle  right)  produces  similar  looking  waveforms,  leaving  relatively 
little  energy  in  the  residual  image  (middle  left).  The  eigenvectors  for  the  larger  eigenvalues  (lower 
left)  and  their  corresponding  reconstructions  (lower  right)  show  the  long-period  noise  character  of 
the  first  4  eigenvectors. 
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Figure  8.  a,b)  Seismic  data  windowed  around  the  P  arrival  for  the  a)  shallow  (2  km)  and  deep  (12 
km)  earthquakes  from  Northridge  (S.  California)  cluster.  The  shallow  event  is  used  as  the 
reference  signal  and  the  deeper  event  is  used  as  the  primary  signal;  noise  has  been  added  to  the 
latter  trace,  c)  Zero-lag  cross-correlation  between  the  primary  and  CANC-filtered  traces  for  a 
range  of  signal  alignments  (see  text),  d)  SNR  for  the  CANC-filtered  trace  versus  the  error  trace  for 
a  range  of  signal  alignments,  e)  Common  signal  obtained  from  CANC  adaptive  filter,  f)  Error 
trace.  Coherency  filtering  applied  to  data  in  a)  and  b)  results  in  the  signal  traces  displayed  in  g) 
and  h),  respectively. 
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Figure  9.  A  horizontal  slice  through 
the  true  synthetic  velocity  model.  The 
true  velocity  model  in  3D  is  similar  to 
a  ’’vertical  sandwich,”  with  velocity 
constant  with  depth. 
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Figure  10.  Event  locations  (filled  circles)  and  stations  (triangles) 
used  for  the  synthetic  data  set.  The  inversion  grid  used  in  the 
standard  and  DD  tomography  solutions  is  shown  as  crosses.  The 
inversion  grid  points  are  at  X=-35,  -15,  0,  2,  4,  6,  20,  35  km,  at 
Y=-60,  -40  -20,  0,  20,  40  km,  and  at  Z=0,  3,  7,  1 1,  16  km. 


Figure  1 1.  Horizontal  slices  through  the  velocity  models  trom  (a)  standard  tomography  and  (b) 
DD  tomography,  and  the  difference  between  (c)  the  DD  tomography  solution  and  the  true  model 
and  (d)  the  standard  tomography  solution  and  the  true  model.  Black  dots  indicate  the  earthquake 
hypocenters  within  half  the  grid-size  of  the  slice. 
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Figure  12.  Comparison  of  earthquake  locations  along  the  Hayward  fault  determined  using  (a)  DD 
tomography,  (b)  DD  location,  (c)  standard  tomography,  and  (d)  catalog  data  (no  DD  data).  Each 
part  shows  a  lat-lon  plot  (upper  left)  and  zoomed-in  plots  (rotated  to  fault-parallel/fault-normal) 
of  epicenters  (upper  right)  and  fault-normal  (lower  left)  and  fault-parallel  (lower  right)  depth 
sections.  Note  how  closely  the  hypocenters  in  panels  a  and  b  resemble  each  other,  but  with  a 
systematic  SW  shift.  The  tight  clustering  is  not  evident  in  either  the  standard  tomography 
solution  (c)  or  the  catalog  locations  using  absolute  picks  only  (d). 
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Figure  13.  (a)  The  left-hand  panels  show  fault-normal  slices  through  the  DD  tomography  model. 
The  model  has  sharper  features  and  a  better  correspondence  between  the  main  region  of 
seismicity  and  the  position  of  the  sharp  velocity  contrast  compared  to  the  standard  results  in  b. 
(b)  The  right-hand  panels  show  fault-normal  slices  through  the  standard  tomography  model.  The 
model  has  features  that  are  less  sharp  and  a  poorer  correspondence  between  the  main  region  of 
seismicity  and  the  position  of  the  velocity  contrast  than  the  DD  results 


29 


Table  1.  The  absolute  differences  between  the  true  locations  and  those  from  the  DD  location 
method  based  on  1 D  velocity  model,  standard  tomography,  and  DD  tomography. 


Median  value  (km) 

Standard  deviation  (km) 

Latitude 

Longitude 

Depth 

Latitude 

Longitude 

Depth 

DD  location  ( 1 D) 

1.131 

1.235 

1.123 

0.976 

0.941 

1.658 

Standard  tomography 

0.320 

0.295 

0.460 

0.399 

0.342 

0.575 

DD  tomography 

0.238 

0.218 

0.329 

0.288 

0.314 

0.427 
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List  of  Symbols,  Abbreviations,  and  Acronyms 

relative  time  delay  in  lag  number  for  an  event  pair 


A  k(sub) 

^ij 

A  k(bsl) 

^ij 

A  k(bs2) 

^ij 

ID 

3D 

BCSEIS 

BS 

CANC 

CC 

cc'’* 

CCij'' 

CCij'"“ 

ccsub 

cc'im' 

Cc'ini2 

Cc'im' 

CSS-3.0 

DD 

ds 

GT 

hypoDD 

'J 

k 

KLD 

LI 

M 

NCEDC 

NMT 

PLRR 

RMS 

SNR 

T 

relative  time  delay  in  fraction  of  a  sample  obtained  by  the  cross  spectrum  phase 
bispectrum  time  delay  in  lag  number  using  the  band-pass  filtered  waveforms 
bispectrum  time  delay  in  lag  number  using  the  raw  waveforms 
one-dimensional 
three-dimensional 

software  package  for  bispectrum  verification  of  cross-correlation  data 
bispectrum 

Correlated  data  Adaptive  Noise  Canceling 
cross-correlation 

threshold  for  performing  the  bispectrum  time  delay  estimation 

cross-correlation  value  for  event  pair  i  and  j  at  station  k 

maximum  cross-correlation  value  for  an  event  pair 

threshold  for  performing  sub-sample  time  delay  calculation 

threshold  value  to  select  cross-correlation  time  delay  estimates 

threshold  value  to  carry  out  bispectrum  verification  for  an  event  pair 

threshold  value  to  carry  out  bispectrum  verification  for  an  individual  station 

Center  for  Seismic  Studies  format  version  3.0 

double-difference 

an  element  of  path  length 

ground  truth 

double-difference  location  algorithm 
event  numbers 
station  number 

Karhunen-Loeve  decomposition 
one-norm  (absolute  value)  measure  of  size 
magnitude 

Northern  California  Earthquake  Data  Center 

New  Mexico  Tech 

Pick-Locate-Repick-Relocate 
root-mean-square 
signal  to  noise  ratio 
body  wave  arrival  time 

x' 

tomoDD 

u 

UC 

UW 

WCC 

XI,  X2,  X3 

X,  Y,Z 

origin  time  of  event  i 
double-difference  tomography  algorithm 
slowness  field  (reciprocal  of  velocity) 

University  of  California 

University  of  Wisconsin-Madison 
waveform  cross-correlation 
source  coordinates 

coordinate  axes 
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